!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine ucs (ro, mx, my, mz, bx, by, bz, en, gx, gy,                                                                            &
   nx, ny, ng, n_neighs, dx, dy, dt, gm,                                                                                           &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd,                                                 &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny, ng, n_neighs
   real(8), intent(in) :: dx, dy, dt, gm
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy
   real(8), dimension(ny, n_neighs), optional, intent(inout) ::                                                                    &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd
   real(8), dimension(nx, n_neighs), optional, intent(inout) ::                                                                    &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd

   integer(4) :: i, ip1, im1, j, jp1, jm1

   real(8), dimension(nx, ny) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre
   real(8), dimension(nx, ny) :: ro_l_x, mx_l_x, my_l_x, mz_l_x, bx_l_x, by_l_x, bz_l_x, en_l_x
   real(8), dimension(nx, ny) :: ro_l_y, mx_l_y, my_l_y, mz_l_y, bx_l_y, by_l_y, bz_l_y, en_l_y
   real(8), dimension(nx, ny) :: ro_r_x, mx_r_x, my_r_x, mz_r_x, bx_r_x, by_r_x, bz_r_x, en_r_x
   real(8), dimension(nx, ny) :: ro_r_y, mx_r_y, my_r_y, mz_r_y, bx_r_y, by_r_y, bz_r_y, en_r_y
   real(8), dimension(nx, ny) :: ro_mod_x, mx_mod_x, my_mod_x, mz_mod_x, bx_mod_x, by_mod_x, bz_mod_x, en_mod_x
   real(8), dimension(nx, ny) :: ro_mod_y, mx_mod_y, my_mod_y, mz_mod_y, bx_mod_y, by_mod_y, bz_mod_y, en_mod_y
   real(8), dimension(nx, ny) :: ro_g1, mx_g1, my_g1, mz_g1, bx_g1, by_g1, bz_g1, en_g1
   real(8), dimension(nx, ny) :: fxro_pre, fxmx_pre, fxmy_pre, fxmz_pre, fxbx_pre, fxby_pre, fxbz_pre, fxen_pre
   real(8), dimension(nx, ny) :: fyro_pre, fymx_pre, fymy_pre, fymz_pre, fybx_pre, fyby_pre, fybz_pre, fyen_pre
   real(8), dimension(nx, ny) :: emf, emf_g1

   real(8) :: fxro_l, fxmx_l, fxmy_l, fxmz_l, fxbx_l, fxby_l, fxbz_l, fxen_l
   real(8) :: fxro_r, fxmx_r, fxmy_r, fxmz_r, fxbx_r, fxby_r, fxbz_r, fxen_r
   real(8) :: fyro_l, fymx_l, fymy_l, fymz_l, fybx_l, fyby_l, fybz_l, fyen_l
   real(8) :: fyro_r, fymx_r, fymy_r, fymz_r, fybx_r, fyby_r, fybz_r, fyen_r
   real(8) :: limiter
   integer(4) :: sgn
   real(8) :: dt_dx, dt_dy, u1, u2, u3, vx, vy, vz, b2, v2, pr, gmm1, theta

!----------------------------------------------------------------------------------------------------------------------------------|
   dt_dx = dt / dx
   dt_dy = dt / dy
   gmm1 = gm - 1.0d0
   theta = 1.0d0

   do j = 2, ny - 1
   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1
      jp1 = j + 1
      jm1 = j - 1

      u1 = (ro(i, j) - ro(im1, j)) * theta
      u2 = (ro(ip1, j) - ro(im1, j)) * 0.5d0
      u3 = (ro(ip1, j) - ro(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_l_x(i, j) = ro(i, j) + 0.5d0 * limiter
      ro_r_x(i, j) = ro(i, j) - 0.5d0 * limiter
      ro_mod_x(i, j) = limiter

      u1 = (ro(i, j) - ro(i, jm1)) * theta
      u2 = (ro(i, jp1) - ro(i, jm1)) * 0.5d0
      u3 = (ro(i, jp1) - ro(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_l_y(i, j) = ro(i, j) + 0.5d0 * limiter
      ro_r_y(i, j) = ro(i, j) - 0.5d0 * limiter
      ro_mod_y(i, j) = limiter

      u1 = (mx(i, j) - mx(im1, j)) * theta
      u2 = (mx(ip1, j) - mx(im1, j)) * 0.5d0
      u3 = (mx(ip1, j) - mx(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_l_x(i, j) = mx(i, j) + 0.5d0 * limiter
      mx_r_x(i, j) = mx(i, j) - 0.5d0 * limiter
      mx_mod_x(i, j) = limiter

      u1 = (mx(i, j) - mx(i, jm1)) * theta
      u2 = (mx(i, jp1) - mx(i, jm1)) * 0.5d0
      u3 = (mx(i, jp1) - mx(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_l_y(i, j) = mx(i, j) + 0.5d0 * limiter
      mx_r_y(i, j) = mx(i, j) - 0.5d0 * limiter
      mx_mod_y(i, j) = limiter

      u1 = (my(i, j) - my(im1, j)) * theta
      u2 = (my(ip1, j) - my(im1, j)) * 0.5d0
      u3 = (my(ip1, j) - my(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_l_x(i, j) = my(i, j) + 0.5d0 * limiter
      my_r_x(i, j) = my(i, j) - 0.5d0 * limiter
      my_mod_x(i, j) = limiter

      u1 = (my(i, j) - my(i, jm1)) * theta
      u2 = (my(i, jp1) - my(i, jm1)) * 0.5d0
      u3 = (my(i, jp1) - my(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_l_y(i, j) = my(i, j) + 0.5d0 * limiter
      my_r_y(i, j) = my(i, j) - 0.5d0 * limiter
      my_mod_y(i, j) = limiter

      u1 = (mz(i, j) - mz(im1, j)) * theta
      u2 = (mz(ip1, j) - mz(im1, j)) * 0.5d0
      u3 = (mz(ip1, j) - mz(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_l_x(i, j) = mz(i, j) + 0.5d0 * limiter
      mz_r_x(i, j) = mz(i, j) - 0.5d0 * limiter
      mz_mod_x(i, j) = limiter

      u1 = (mz(i, j) - mz(i, jm1)) * theta
      u2 = (mz(i, jp1) - mz(i, jm1)) * 0.5d0
      u3 = (mz(i, jp1) - mz(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_l_y(i, j) = mz(i, j) + 0.5d0 * limiter
      mz_r_y(i, j) = mz(i, j) - 0.5d0 * limiter
      mz_mod_y(i, j) = limiter

      u1 = (bx(i, j) - bx(im1, j)) * theta
      u2 = (bx(ip1, j) - bx(im1, j)) * 0.5d0
      u3 = (bx(ip1, j) - bx(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_l_x(i, j) = bx(i, j) + 0.5d0 * limiter
      bx_r_x(i, j) = bx(i, j) - 0.5d0 * limiter
      bx_mod_x(i, j) = limiter

      u1 = (bx(i, j) - bx(i, jm1)) * theta
      u2 = (bx(i, jp1) - bx(i, jm1)) * 0.5d0
      u3 = (bx(i, jp1) - bx(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_l_y(i, j) = bx(i, j) + 0.5d0 * limiter
      bx_r_y(i, j) = bx(i, j) - 0.5d0 * limiter
      bx_mod_y(i, j) = limiter

      u1 = (by(i, j) - by(im1, j)) * theta
      u2 = (by(ip1, j) - by(im1, j)) * 0.5d0
      u3 = (by(ip1, j) - by(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_l_x(i, j) = by(i, j) + 0.5d0 * limiter
      by_r_x(i, j) = by(i, j) - 0.5d0 * limiter
      by_mod_x(i, j) = limiter

      u1 = (by(i, j) - by(i, jm1)) * theta
      u2 = (by(i, jp1) - by(i, jm1)) * 0.5d0
      u3 = (by(i, jp1) - by(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_l_y(i, j) = by(i, j) + 0.5d0 * limiter
      by_r_y(i, j) = by(i, j) - 0.5d0 * limiter
      by_mod_y(i, j) = limiter

      u1 = (bz(i, j) - bz(im1, j)) * theta
      u2 = (bz(ip1, j) - bz(im1, j)) * 0.5d0
      u3 = (bz(ip1, j) - bz(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_l_x(i, j) = bz(i, j) + 0.5d0 * limiter
      bz_r_x(i, j) = bz(i, j) - 0.5d0 * limiter
      bz_mod_x(i, j) = limiter

      u1 = (bz(i, j) - bz(i, jm1)) * theta
      u2 = (bz(i, jp1) - bz(i, jm1)) * 0.5d0
      u3 = (bz(i, jp1) - bz(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_l_y(i, j) = bz(i, j) + 0.5d0 * limiter
      bz_r_y(i, j) = bz(i, j) - 0.5d0 * limiter
      bz_mod_y(i, j) = limiter

      u1 = (en(i, j) - en(im1, j)) * theta
      u2 = (en(ip1, j) - en(im1, j)) * 0.5d0
      u3 = (en(ip1, j) - en(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_l_x(i, j) = en(i, j) + 0.5d0 * limiter
      en_r_x(i, j) = en(i, j) - 0.5d0 * limiter
      en_mod_x(i, j) = limiter

      u1 = (en(i, j) - en(i, jm1)) * theta
      u2 = (en(i, jp1) - en(i, jm1)) * 0.5d0
      u3 = (en(i, jp1) - en(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_l_y(i, j) = en(i, j) + 0.5d0 * limiter
      en_r_y(i, j) = en(i, j) - 0.5d0 * limiter
      en_mod_y(i, j) = limiter

   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      vx = mx_l_x(i, j) / ro_l_x(i, j)
      vy = my_l_x(i, j) / ro_l_x(i, j)
      vz = mz_l_x(i, j) / ro_l_x(i, j)
      b2 = bx_l_x(i, j) * bx_l_x(i, j) + by_l_x(i, j) * by_l_x(i, j) + bz_l_x(i, j) * bz_l_x(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l_x(i, j) - v2 * ro_l_x(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      fxro_l = mx_l_x(i, j)
	   fxmx_l = vx * mx_l_x(i, j) + pr + b2 / 2.0d0 - bx_l_x(i, j) * bx_l_x(i, j)
	   fxmy_l = vx * my_l_x(i, j) - bx_l_x(i, j) * by_l_x(i, j)
	   fxmz_l = vx * mz_l_x(i, j) - bx_l_x(i, j) * bz_l_x(i, j)
	   fxbx_l = 0.0d0
	   fxby_l = vx * by_l_x(i, j) - bx_l_x(i, j) * vy
	   fxbz_l = vx * bz_l_x(i, j) - bx_l_x(i, j) * vz
	   fxen_l = (en_l_x(i, j) + pr + b2 / 2.0d0) * vx - (bx_l_x(i, j) * vx + by_l_x(i, j) * vy + bz_l_x(i, j) * vz) * bx_l_x(i, j)
      vx = mx_r_x(i, j) / ro_r_x(i, j)
      vy = my_r_x(i, j) / ro_r_x(i, j)
      vz = mz_r_x(i, j) / ro_r_x(i, j)
      b2 = bx_r_x(i, j) * bx_r_x(i, j) + by_r_x(i, j) * by_r_x(i, j) + bz_r_x(i, j) * bz_r_x(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r_x(i, j) - v2 * ro_r_x(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      fxro_r = mx_r_x(i, j)
	   fxmx_r = vx * mx_r_x(i, j) + pr + b2 / 2.0d0 - bx_r_x(i, j) * bx_r_x(i, j)
	   fxmy_r = vx * my_r_x(i, j) - bx_r_x(i, j) * by_r_x(i, j)
	   fxmz_r = vx * mz_r_x(i, j) - bx_r_x(i, j) * bz_r_x(i, j)
	   fxbx_r = 0.0d0
	   fxby_r = vx * by_r_x(i, j) - bx_r_x(i, j) * vy
	   fxbz_r = vx * bz_r_x(i, j) - bx_r_x(i, j) * vz
	   fxen_r = (en_r_x(i, j) + pr + b2 / 2.0d0) * vx - (bx_r_x(i, j) * vx + by_r_x(i, j) * vy + bz_r_x(i, j) * vz) * bx_r_x(i, j)

      ro_pre(i, j) = ro(i, j) - 0.5d0 * dt_dx * (fxro_l - fxro_r)
      mx_pre(i, j) = mx(i, j) - 0.5d0 * dt_dx * (fxmx_l - fxmx_r)
      my_pre(i, j) = my(i, j) - 0.5d0 * dt_dx * (fxmy_l - fxmy_r)
      mz_pre(i, j) = mz(i, j) - 0.5d0 * dt_dx * (fxmz_l - fxmz_r)
      bx_pre(i, j) = bx(i, j) - 0.5d0 * dt_dx * (fxbx_l - fxbx_r)
      by_pre(i, j) = by(i, j) - 0.5d0 * dt_dx * (fxby_l - fxby_r)
      bz_pre(i, j) = bz(i, j) - 0.5d0 * dt_dx * (fxbz_l - fxbz_r)
      en_pre(i, j) = en(i, j) - 0.5d0 * dt_dx * (fxen_l - fxen_r)

      vx = mx_l_y(i, j) / ro_l_y(i, j)
      vy = my_l_y(i, j) / ro_l_y(i, j)
      vz = mz_l_y(i, j) / ro_l_y(i, j)
      b2 = bx_l_y(i, j) * bx_l_y(i, j) + by_l_y(i, j) * by_l_y(i, j) + bz_l_y(i, j) * bz_l_y(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l_y(i, j) - v2 * ro_l_y(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      fyro_l = my_l_y(i, j)
	   fymx_l = vy * mx_l_y(i, j) - by_l_y(i, j) * bx_l_y(i, j)
	   fymy_l = vy * my_l_y(i, j) + pr + b2 / 2.0d0 - by_l_y(i, j) * by_l_y(i, j)
	   fymz_l = vy * mz_l_y(i, j) - by_l_y(i, j) * bz_l_y(i, j)
	   fybx_l = vy * bx_l_y(i, j) - by_l_y(i, j) * vx
	   fyby_l = 0.0d0
	   fybz_l = vy * bz_l_y(i, j) - by_l_y(i, j) * vz
	   fyen_l = (en_l_y(i, j) + pr + b2 / 2.0d0) * vy - (bx_l_y(i, j) * vx + by_l_y(i, j) * vy + bz_l_y(i, j) * vz) * by_l_y(i, j)
      vx = mx_r_y(i, j) / ro_r_y(i, j)
      vy = my_r_y(i, j) / ro_r_y(i, j)
      vz = mz_r_y(i, j) / ro_r_y(i, j)
      b2 = bx_r_y(i, j) * bx_r_y(i, j) + by_r_y(i, j) * by_r_y(i, j) + bz_r_y(i, j) * bz_r_y(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r_y(i, j) - v2 * ro_r_y(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      fyro_r = my_r_y(i, j)
	   fymx_r = vy * mx_r_y(i, j) - by_r_y(i, j) * bx_r_y(i, j)
	   fymy_r = vy * my_r_y(i, j) + pr + b2 / 2.0d0 - by_r_y(i, j) * by_r_y(i, j)
	   fymz_r = vy * mz_r_y(i, j) - by_r_y(i, j) * bz_r_y(i, j)
	   fybx_r = vy * bx_r_y(i, j) - by_r_y(i, j) * vx
	   fyby_r = 0.0d0
	   fybz_r = vy * bz_r_y(i, j) - by_r_y(i, j) * vz
	   fyen_r = (en_r_y(i, j) + pr + b2 / 2.0d0) * vy - (bx_r_y(i, j) * vx + by_r_y(i, j) * vy + bz_r_y(i, j) * vz) * by_r_y(i, j)

      ro_pre(i, j) = ro_pre(i, j) - 0.5d0 * dt_dy * (fyro_l - fyro_r)
      mx_pre(i, j) = mx_pre(i, j) - 0.5d0 * dt_dy * (fymx_l - fymx_r)
      my_pre(i, j) = my_pre(i, j) - 0.5d0 * dt_dy * (fymy_l - fymy_r)
      mz_pre(i, j) = mz_pre(i, j) - 0.5d0 * dt_dy * (fymz_l - fymz_r)
      bx_pre(i, j) = bx_pre(i, j) - 0.5d0 * dt_dy * (fybx_l - fybx_r)
      by_pre(i, j) = by_pre(i, j) - 0.5d0 * dt_dy * (fyby_l - fyby_r)
      bz_pre(i, j) = bz_pre(i, j) - 0.5d0 * dt_dy * (fybz_l - fybz_r)
      en_pre(i, j) = en_pre(i, j) - 0.5d0 * dt_dy * (fyen_l - fyen_r)

      vx = mx_pre(i, j) / ro_pre(i, j)
      vy = my_pre(i, j) / ro_pre(i, j)
      vz = mz_pre(i, j) / ro_pre(i, j)
      b2 = bx_pre(i, j) * bx_pre(i, j) + by_pre(i, j) * by_pre(i, j) + bz_pre(i, j) * bz_pre(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_pre(i, j) - v2 * ro_pre(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      fxro_pre(i, j) = mx_pre(i, j)
	   fxmx_pre(i, j) = vx * mx_pre(i, j) + pr + b2 / 2.0d0 - bx_pre(i, j) * bx_pre(i, j)
	   fxmy_pre(i, j) = vx * my_pre(i, j) - bx_pre(i, j) * by_pre(i, j)
	   fxmz_pre(i, j) = vx * mz_pre(i, j) - bx_pre(i, j) * bz_pre(i, j)
	   fxbx_pre(i, j) = 0.0d0
	   fxby_pre(i, j) = vx * by_pre(i, j) - bx_pre(i, j) * vy
	   fxbz_pre(i, j) = vx * bz_pre(i, j) - bx_pre(i, j) * vz
	   fxen_pre(i, j) = (en_pre(i, j) + pr + b2 / 2.0d0) * vx - (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * bx_pre(i, j)
      fyro_pre(i, j) = my_pre(i, j)
	   fymx_pre(i, j) = vy * mx_pre(i, j) - by_pre(i, j) * bx_pre(i, j)
	   fymy_pre(i, j) = vy * my_pre(i, j) + pr + b2 / 2.0d0 - by_pre(i, j) * by_pre(i, j)
	   fymz_pre(i, j) = vy * mz_pre(i, j) - by_pre(i, j) * bz_pre(i, j)
	   fybx_pre(i, j) = vy * bx_pre(i, j) - by_pre(i, j) * vx
	   fyby_pre(i, j) = 0.0d0
	   fybz_pre(i, j) = vy * bz_pre(i, j) - by_pre(i, j) * vz
	   fyen_pre(i, j) = (en_pre(i, j) + pr + b2 / 2.0d0) * vy - (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * by_pre(i, j)

      emf(i, j) = (bx(i, j) * my(i, j) - by(i, j) * mx(i, j)) / ro(i, j)
   enddo
   enddo

   do j = 1, ny - 1
   do i = 1, nx - 1
      ip1 = i + 1
      jp1 = j + 1

      ro_g1(i, j) = 0.25d0 * (ro(i, j) + ro(ip1, j) + ro(ip1, jp1) + ro(i, jp1) +                                                  &
         0.25d0 * (ro_mod_x(i, j) - ro_mod_x(ip1, j) + ro_mod_x(i, jp1) - ro_mod_x(ip1, jp1) +                                     &
                   ro_mod_y(i, j) - ro_mod_y(i, jp1) + ro_mod_y(ip1, j) - ro_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxro_pre(ip1, j) - fxro_pre(i, j) + fxro_pre(ip1, jp1) - fxro_pre(i, jp1)) +                            &
                  dt_dy * (fyro_pre(i, jp1) - fyro_pre(i, j) + fyro_pre(ip1, jp1) - fyro_pre(ip1, j)))
      mx_g1(i, j) = 0.25d0 * (mx(i, j) + mx(ip1, j) + mx(ip1, jp1) + mx(i, jp1) +                                                  &
         0.25d0 * (mx_mod_x(i, j) - mx_mod_x(ip1, j) + mx_mod_x(i, jp1) - mx_mod_x(ip1, jp1) +                                     &
                   mx_mod_y(i, j) - mx_mod_y(i, jp1) + mx_mod_y(ip1, j) - mx_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxmx_pre(ip1, j) - fxmx_pre(i, j) + fxmx_pre(ip1, jp1) - fxmx_pre(i, jp1)) +                            &
                  dt_dy * (fymx_pre(i, jp1) - fymx_pre(i, j) + fymx_pre(ip1, jp1) - fymx_pre(ip1, j))) +                           &
         dt * 0.25d0 * (ro(i, j) * gx(i, j) + ro(ip1, j) * gx(ip1, j) + ro(ip1, jp1) * gx(ip1, jp1) + ro(i, jp1) * gx(i, jp1))
      my_g1(i, j) = 0.25d0 * (my(i, j) + my(ip1, j) + my(ip1, jp1) + my(i, jp1) +                                                  &
         0.25d0 * (my_mod_x(i, j) - my_mod_x(ip1, j) + my_mod_x(i, jp1) - my_mod_x(ip1, jp1) +                                     &
                   my_mod_y(i, j) - my_mod_y(i, jp1) + my_mod_y(ip1, j) - my_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxmy_pre(ip1, j) - fxmy_pre(i, j) + fxmy_pre(ip1, jp1) - fxmy_pre(i, jp1)) +                            &
                  dt_dy * (fymy_pre(i, jp1) - fymy_pre(i, j) + fymy_pre(ip1, jp1) - fymy_pre(ip1, j))) +                           &
         dt * 0.25d0 * (ro(i, j) * gy(i, j) + ro(ip1, j) * gy(ip1, j) + ro(ip1, jp1) * gy(ip1, jp1) + ro(i, jp1) * gy(i, jp1))
      mz_g1(i, j) = 0.25d0 * (mz(i, j) + mz(ip1, j) + mz(ip1, jp1) + mz(i, jp1) +                                                  &
         0.25d0 * (mz_mod_x(i, j) - mz_mod_x(ip1, j) + mz_mod_x(i, jp1) - mz_mod_x(ip1, jp1) +                                     &
                   mz_mod_y(i, j) - mz_mod_y(i, jp1) + mz_mod_y(ip1, j) - mz_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxmz_pre(ip1, j) - fxmz_pre(i, j) + fxmz_pre(ip1, jp1) - fxmz_pre(i, jp1)) +                            &
                  dt_dy * (fymz_pre(i, jp1) - fymz_pre(i, j) + fymz_pre(ip1, jp1) - fymz_pre(ip1, j)))
      bx_g1(i, j) = 0.25d0 * (bx(i, j) + bx(ip1, j) + bx(ip1, jp1) + bx(i, jp1) +                                                  &
         0.25d0 * (bx_mod_x(i, j) - bx_mod_x(ip1, j) + bx_mod_x(i, jp1) - bx_mod_x(ip1, jp1) +                                     &
                   bx_mod_y(i, j) - bx_mod_y(i, jp1) + bx_mod_y(ip1, j) - bx_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxbx_pre(ip1, j) - fxbx_pre(i, j) + fxbx_pre(ip1, jp1) - fxbx_pre(i, jp1)) +                            &
                  dt_dy * (fybx_pre(i, jp1) - fybx_pre(i, j) + fybx_pre(ip1, jp1) - fybx_pre(ip1, j)))
      by_g1(i, j) = 0.25d0 * (by(i, j) + by(ip1, j) + by(ip1, jp1) + by(i, jp1) +                                                  &
         0.25d0 * (by_mod_x(i, j) - by_mod_x(ip1, j) + by_mod_x(i, jp1) - by_mod_x(ip1, jp1) +                                     &
                   by_mod_y(i, j) - by_mod_y(i, jp1) + by_mod_y(ip1, j) - by_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxby_pre(ip1, j) - fxby_pre(i, j) + fxby_pre(ip1, jp1) - fxby_pre(i, jp1)) +                            &
                  dt_dy * (fyby_pre(i, jp1) - fyby_pre(i, j) + fyby_pre(ip1, jp1) - fyby_pre(ip1, j)))
      bz_g1(i, j) = 0.25d0 * (bz(i, j) + bz(ip1, j) + bz(ip1, jp1) + bz(i, jp1) +                                                  &
         0.25d0 * (bz_mod_x(i, j) - bz_mod_x(ip1, j) + bz_mod_x(i, jp1) - bz_mod_x(ip1, jp1) +                                     &
                   bz_mod_y(i, j) - bz_mod_y(i, jp1) + bz_mod_y(ip1, j) - bz_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxbz_pre(ip1, j) - fxbz_pre(i, j) + fxbz_pre(ip1, jp1) - fxbz_pre(i, jp1)) +                            &
                  dt_dy * (fybz_pre(i, jp1) - fybz_pre(i, j) + fybz_pre(ip1, jp1) - fybz_pre(ip1, j)))
      en_g1(i, j) = 0.25d0 * (en(i, j) + en(ip1, j) + en(ip1, jp1) + en(i, jp1) +                                                  &
         0.25d0 * (en_mod_x(i, j) - en_mod_x(ip1, j) + en_mod_x(i, jp1) - en_mod_x(ip1, jp1) +                                     &
                   en_mod_y(i, j) - en_mod_y(i, jp1) + en_mod_y(ip1, j) - en_mod_y(ip1, jp1))) -                                   &
         0.5d0 * (dt_dx * (fxen_pre(ip1, j) - fxen_pre(i, j) + fxen_pre(ip1, jp1) - fxen_pre(i, jp1)) +                            &
                  dt_dy * (fyen_pre(i, jp1) - fyen_pre(i, j) + fyen_pre(ip1, jp1) - fyen_pre(ip1, j))) +                           &
         dt * 0.25d0 * (mx(i, j) * gx(i, j) + mx(ip1, j) * gx(ip1, j) + mx(ip1, jp1) * gx(ip1, jp1) + mx(i, jp1) * gx(i, jp1)) +   &
         dt * 0.25d0 * (my(i, j) * gy(i, j) + my(ip1, j) * gy(ip1, j) + my(ip1, jp1) * gy(ip1, jp1) + my(i, jp1) * gy(i, jp1))

      emf_g1(i, j) = (bx_g1(i, j) * my_g1(i, j) - by_g1(i, j) * mx_g1(i, j)) / ro_g1(i, j)
      emf_g1(i, j) = 0.5d0 * (emf_g1(i, j) + 0.25d0 * (emf(i, j) + emf(i, jp1) + emf(ip1, j) + emf(ip1, jp1)))
   enddo
   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1
      jp1 = j + 1
      jm1 = j - 1

      u1 = (ro_g1(i, j) - ro_g1(im1, j)) * theta
      u2 = (ro_g1(ip1, j) - ro_g1(im1, j)) * 0.5d0
      u3 = (ro_g1(ip1, j) - ro_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_mod_x(i, j) = limiter

      u1 = (ro_g1(i, j) - ro_g1(i, jm1)) * theta
      u2 = (ro_g1(i, jp1) - ro_g1(i, jm1)) * 0.5d0
      u3 = (ro_g1(i, jp1) - ro_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_mod_y(i, j) = limiter

      u1 = (mx_g1(i, j) - mx_g1(im1, j)) * theta
      u2 = (mx_g1(ip1, j) - mx_g1(im1, j)) * 0.5d0
      u3 = (mx_g1(ip1, j) - mx_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_mod_x(i, j) = limiter

      u1 = (mx_g1(i, j) - mx_g1(i, jm1)) * theta
      u2 = (mx_g1(i, jp1) - mx_g1(i, jm1)) * 0.5d0
      u3 = (mx_g1(i, jp1) - mx_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_mod_y(i, j) = limiter

      u1 = (my_g1(i, j) - my_g1(im1, j)) * theta
      u2 = (my_g1(ip1, j) - my_g1(im1, j)) * 0.5d0
      u3 = (my_g1(ip1, j) - my_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_mod_x(i, j) = limiter

      u1 = (my_g1(i, j) - my_g1(i, jm1)) * theta
      u2 = (my_g1(i, jp1) - my_g1(i, jm1)) * 0.5d0
      u3 = (my_g1(i, jp1) - my_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_mod_y(i, j) = limiter

      u1 = (mz_g1(i, j) - mz_g1(im1, j)) * theta
      u2 = (mz_g1(ip1, j) - mz_g1(im1, j)) * 0.5d0
      u3 = (mz_g1(ip1, j) - mz_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_mod_x(i, j) = limiter

      u1 = (mz_g1(i, j) - mz_g1(i, jm1)) * theta
      u2 = (mz_g1(i, jp1) - mz_g1(i, jm1)) * 0.5d0
      u3 = (mz_g1(i, jp1) - mz_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_mod_y(i, j) = limiter

      u1 = (bx_g1(i, j) - bx_g1(im1, j)) * theta
      u2 = (bx_g1(ip1, j) - bx_g1(im1, j)) * 0.5d0
      u3 = (bx_g1(ip1, j) - bx_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_mod_x(i, j) = limiter

      u1 = (bx_g1(i, j) - bx_g1(i, jm1)) * theta
      u2 = (bx_g1(i, jp1) - bx_g1(i, jm1)) * 0.5d0
      u3 = (bx_g1(i, jp1) - bx_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_mod_y(i, j) = limiter

      u1 = (by_g1(i, j) - by_g1(im1, j)) * theta
      u2 = (by_g1(ip1, j) - by_g1(im1, j)) * 0.5d0
      u3 = (by_g1(ip1, j) - by_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_mod_x(i, j) = limiter

      u1 = (by_g1(i, j) - by_g1(i, jm1)) * theta
      u2 = (by_g1(i, jp1) - by_g1(i, jm1)) * 0.5d0
      u3 = (by_g1(i, jp1) - by_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_mod_y(i, j) = limiter

      u1 = (bz_g1(i, j) - bz_g1(im1, j)) * theta
      u2 = (bz_g1(ip1, j) - bz_g1(im1, j)) * 0.5d0
      u3 = (bz_g1(ip1, j) - bz_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_mod_x(i, j) = limiter

      u1 = (bz_g1(i, j) - bz_g1(i, jm1)) * theta
      u2 = (bz_g1(i, jp1) - bz_g1(i, jm1)) * 0.5d0
      u3 = (bz_g1(i, jp1) - bz_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_mod_y(i, j) = limiter

      u1 = (en_g1(i, j) - en_g1(im1, j)) * theta
      u2 = (en_g1(ip1, j) - en_g1(im1, j)) * 0.5d0
      u3 = (en_g1(ip1, j) - en_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_mod_x(i, j) = limiter

      u1 = (en_g1(i, j) - en_g1(i, jm1)) * theta
      u2 = (en_g1(i, jp1) - en_g1(i, jm1)) * 0.5d0
      u3 = (en_g1(i, jp1) - en_g1(i, j)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_mod_y(i, j) = limiter

!      bx_g1(i, j) = 0.25d0 * (bx(i, j) + bx(ip1, j) + bx(ip1, jp1) + bx(i, jp1)) - 0.5d0 * dt_dy * (emf_g1(i, jp1) - emf_g1(i, jm1))
!      by_g1(i, j) = 0.25d0 * (by(i, j) + by(ip1, j) + by(ip1, jp1) + by(i, jp1)) + 0.5d0 * dt_dx * (emf_g1(ip1, j) - emf_g1(im1, j))
   enddo
   enddo

   do j = 2, ny
   do i = 2, nx
      im1 = i - 1
      jm1 = j - 1
      ro(i, j) = 0.25d0 * (ro_g1(im1, j) + ro_g1(i, j) + ro_g1(im1, jm1) + ro_g1(i, jm1)) +                                        &
         0.0625d0 * (ro_mod_x(im1, j) - ro_mod_x(i, j) + ro_mod_x(im1, jm1) - ro_mod_x(i, jm1) +                                   &
                     ro_mod_y(i, jm1) - ro_mod_y(i, j) + ro_mod_y(im1, jm1) - ro_mod_y(im1, j))
      mx(i, j) = 0.25d0 * (mx_g1(im1, j) + mx_g1(i, j) + mx_g1(im1, jm1) + mx_g1(i, jm1)) +                                        &
         0.0625d0 * (mx_mod_x(im1, j) - mx_mod_x(i, j) + mx_mod_x(im1, jm1) - mx_mod_x(i, jm1) +                                   &
                     mx_mod_y(i, jm1) - mx_mod_y(i, j) + mx_mod_y(im1, jm1) - mx_mod_y(im1, j))
      my(i, j) = 0.25d0 * (my_g1(im1, j) + my_g1(i, j) + my_g1(im1, jm1) + my_g1(i, jm1)) +                                        &
         0.0625d0 * (my_mod_x(im1, j) - my_mod_x(i, j) + my_mod_x(im1, jm1) - my_mod_x(i, jm1) +                                   &
                     my_mod_y(i, jm1) - my_mod_y(i, j) + my_mod_y(im1, jm1) - my_mod_y(im1, j))
      mz(i, j) = 0.25d0 * (mz_g1(im1, j) + mz_g1(i, j) + mz_g1(im1, jm1) + mz_g1(i, jm1)) +                                        &
         0.0625d0 * (mz_mod_x(im1, j) - mz_mod_x(i, j) + mz_mod_x(im1, jm1) - mz_mod_x(i, jm1) +                                   &
                     mz_mod_y(i, jm1) - mz_mod_y(i, j) + mz_mod_y(im1, jm1) - mz_mod_y(im1, j))
      bx(i, j) = 0.25d0 * (bx_g1(im1, j) + bx_g1(i, j) + bx_g1(im1, jm1) + bx_g1(i, jm1)) +                                        &
         0.0625d0 * (bx_mod_x(im1, j) - bx_mod_x(i, j) + bx_mod_x(im1, jm1) - bx_mod_x(i, jm1) +                                   &
                     bx_mod_y(i, jm1) - bx_mod_y(i, j) + bx_mod_y(im1, jm1) - bx_mod_y(im1, j))
      by(i, j) = 0.25d0 * (by_g1(im1, j) + by_g1(i, j) + by_g1(im1, jm1) + by_g1(i, jm1)) +                                        &
         0.0625d0 * (by_mod_x(im1, j) - by_mod_x(i, j) + by_mod_x(im1, jm1) - by_mod_x(i, jm1) +                                   &
                     by_mod_y(i, jm1) - by_mod_y(i, j) + by_mod_y(im1, jm1) - by_mod_y(im1, j))
      bz(i, j) = 0.25d0 * (bz_g1(im1, j) + bz_g1(i, j) + bz_g1(im1, jm1) + bz_g1(i, jm1)) +                                        &
         0.0625d0 * (bz_mod_x(im1, j) - bz_mod_x(i, j) + bz_mod_x(im1, jm1) - bz_mod_x(i, jm1) +                                   &
                     bz_mod_y(i, jm1) - bz_mod_y(i, j) + bz_mod_y(im1, jm1) - bz_mod_y(im1, j))
      en(i, j) = 0.25d0 * (en_g1(im1, j) + en_g1(i, j) + en_g1(im1, jm1) + en_g1(i, jm1)) +                                        &
         0.0625d0 * (en_mod_x(im1, j) - en_mod_x(i, j) + en_mod_x(im1, jm1) - en_mod_x(i, jm1) +                                   &
                     en_mod_y(i, jm1) - en_mod_y(i, j) + en_mod_y(im1, jm1) - en_mod_y(im1, j))

!      bx(i, j) = 0.25d0 * (bx_g1(i, j) + bx_g1(im1, j) + bx_g1(im1, jm1) + bx_g1(i, jm1))
!      by(i, j) = 0.25d0 * (by_g1(i, j) + by_g1(im1, j) + by_g1(im1, jm1) + by_g1(i, jm1))
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value
!----------------------------------------------------------------------------------------------------------------------------------|
   if (present (fxro_bnd)) then
      do j = 1, ny
         fxro_bnd(j, 1) = fxro_pre(ng, j)
         fxro_bnd(j, 2) = fxro_pre(nx - ng, j)
      enddo
   endif

   if (present (fxmx_bnd)) then
      do j = 1, ny
         fxmx_bnd(j, 1) = fxmx_pre(ng, j)
         fxmx_bnd(j, 2) = fxmx_pre(nx - ng, j)
      enddo
   endif

   if (present (fxmy_bnd)) then
      do j = 1, ny
         fxmy_bnd(j, 1) = fxmy_pre(ng, j)
         fxmy_bnd(j, 2) = fxmy_pre(nx - ng, j)
      enddo
   endif

   if (present (fxmz_bnd)) then
      do j = 1, ny
         fxmz_bnd(j, 1) = fxmz_pre(ng, j)
         fxmz_bnd(j, 2) = fxmz_pre(nx - ng, j)
      enddo
   endif

   if (present (fxbx_bnd)) then
      do j = 1, ny
         fxbx_bnd(j, 1) = fxbx_pre(ng, j)
         fxbx_bnd(j, 2) = fxbx_pre(nx - ng, j)
      enddo
   endif

   if (present (fxby_bnd)) then
      do j = 1, ny
         fxby_bnd(j, 1) = fxby_pre(ng, j)
         fxby_bnd(j, 2) = fxby_pre(nx - ng, j)
      enddo
   endif

   if (present (fxbz_bnd)) then
      do j = 1, ny
         fxbz_bnd(j, 1) = fxbz_pre(ng, j)
         fxbz_bnd(j, 2) = fxbz_pre(nx - ng, j)
      enddo
   endif

   if (present (fxen_bnd)) then
      do j = 1, ny
         fxen_bnd(j, 1) = fxen_pre(ng, j)
         fxen_bnd(j, 2) = fxen_pre(nx - ng, j)
      enddo
   endif

   if (present (fyro_bnd)) then
      do i = 1, nx
         fyro_bnd(i, 1) = fyro_pre(i, ng)
         fyro_bnd(i, 2) = fyro_pre(i, ny - ng)
      enddo
   endif

   if (present (fymx_bnd)) then
      do i = 1, nx
         fymx_bnd(i, 1) = fymx_pre(i, ng)
         fymx_bnd(i, 2) = fymx_pre(i, ny - ng)
      enddo
   endif

   if (present (fymy_bnd)) then
      do i = 1, nx
         fymy_bnd(i, 1) = fymy_pre(i, ng)
         fymy_bnd(i, 2) = fymy_pre(i, ny - ng)
      enddo
   endif

   if (present (fymz_bnd)) then
      do i = 1, nx
         fymz_bnd(i, 1) = fymz_pre(i, ng)
         fymz_bnd(i, 2) = fymz_pre(i, ny - ng)
      enddo
   endif

   if (present (fybx_bnd)) then
      do i = 1, nx
         fybx_bnd(i, 1) = fybx_pre(i, ng)
         fybx_bnd(i, 2) = fybx_pre(i, ny - ng)
      enddo
   endif

   if (present (fyby_bnd)) then
      do i = 1, nx
         fyby_bnd(i, 1) = fyby_pre(i, ng)
         fyby_bnd(i, 2) = fyby_pre(i, ny - ng)
      enddo
   endif

   if (present (fybz_bnd)) then
      do i = 1, nx
         fybz_bnd(i, 1) = fybz_pre(i, ng)
         fybz_bnd(i, 2) = fybz_pre(i, ny - ng)
      enddo
   endif

   if (present (fyen_bnd)) then
      do i = 1, nx
         fyen_bnd(i, 1) = fyen_pre(i, ng)
         fyen_bnd(i, 2) = fyen_pre(i, ny - ng)
      enddo
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine ucs
